<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of v_gaussmixk</title>
  <meta name="keywords" content="v_gaussmixk">
  <meta name="description" content="V_GAUSSMIXK approximate Kullback-Leibler divergence between two GMMs + derivatives">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt;  <a href="index.html">v_mfiles</a> &gt; v_gaussmixk.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for v_mfiles&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->

<h1>v_gaussmixk
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>V_GAUSSMIXK approximate Kullback-Leibler divergence between two GMMs + derivatives</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>function [d,klfg]=v_gaussmixk(mf,vf,wf,mg,vg,wg) </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment">V_GAUSSMIXK approximate Kullback-Leibler divergence between two GMMs + derivatives

 Inputs: with kf &amp; kg mixtures, p data dimensions

   mf(kf,p)                mixture means for GMM f
   vf(kf,p) or vf(p,p,kf)  variances (diagonal or full) for GMM f
   wf(kf,1)                weights for GMM f - must sum to 1
   mg(kg,p)                mixture means for GMM g [g=f if mg,vg,wg omitted]
   vg(kg,p) or vg(p,p,kg)  variances (diagonal or full) for GMM g
   wg(kg,1)                weights for GMM g - must sum to 1

 Outputs:
   d             the approximate KL divergence D(f||g)
   klfg(kf,kg)   the exact KL divergence between the components of f and g

 The Kullback-Leibler (KL) divergence, D(f||g), between two distributions,
 f(x) and g(x) is also known as the &quot;relative v_entropy&quot; of f relative to g.
 It is defined as &lt;log(f(x)) - log(g(x))&gt; where &lt;y(x)&gt; denotes the
 expectation with respect to f(x), i.e. &lt;y(x)&gt; = Integral(f(x) y(x) dx).
 The KL divergence is always &gt;=0 and equals 0 if and only if f(x)=g(x)
 almost everywhere. % It is not a distance because it is not symmetric
 between f and g and also does not satisfy the triangle inequality.
 It is normally appropriate for f(x) to be the &quot;true&quot; distribution and
 g(x) to be an approximation to it. See [1].

 This routine calculates the &quot;variational approximation&quot; to the KL divergence
 from [2] that is exact when f and g are single component gaussians and that is zero
 if f=g. However, the approximation may be negative if f and g are different.

 Refs:
 [1]    S. Kullback and R. Leibler. On information and sufficiency.
       Annals of Mathematical Statistics, 22 (1): 79�86, 1951. doi: 10.1214/aoms/1177729694.
 [2]    J. R. Hershey and P. A. Olsen.
       Approximating the kullback leibler divergence between gaussian mixture models.
       In Proc. IEEE Intl Conf. Acoustics, Speech and Signal Processing, volume 4,
       pages IV�317�IV�320, Apr. 2007. doi: 10.1109/ICASSP.2007.366913.</pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
<li><a href="v_logsum.html" class="code" title="function y=v_logsum(x,d,k)">v_logsum</a>	V_LOGSUM v_logsum(x,d,k)=log(sum(k.*exp(x),d))</li></ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 <a name="_sub0" href="#_subfunctions" class="code">function [d,klfg]=v_gaussmixk(mf,vf,wf,mg,vg,wg)</a>
0002 <span class="comment">%V_GAUSSMIXK approximate Kullback-Leibler divergence between two GMMs + derivatives</span>
0003 <span class="comment">%</span>
0004 <span class="comment">% Inputs: with kf &amp; kg mixtures, p data dimensions</span>
0005 <span class="comment">%</span>
0006 <span class="comment">%   mf(kf,p)                mixture means for GMM f</span>
0007 <span class="comment">%   vf(kf,p) or vf(p,p,kf)  variances (diagonal or full) for GMM f</span>
0008 <span class="comment">%   wf(kf,1)                weights for GMM f - must sum to 1</span>
0009 <span class="comment">%   mg(kg,p)                mixture means for GMM g [g=f if mg,vg,wg omitted]</span>
0010 <span class="comment">%   vg(kg,p) or vg(p,p,kg)  variances (diagonal or full) for GMM g</span>
0011 <span class="comment">%   wg(kg,1)                weights for GMM g - must sum to 1</span>
0012 <span class="comment">%</span>
0013 <span class="comment">% Outputs:</span>
0014 <span class="comment">%   d             the approximate KL divergence D(f||g)</span>
0015 <span class="comment">%   klfg(kf,kg)   the exact KL divergence between the components of f and g</span>
0016 <span class="comment">%</span>
0017 <span class="comment">% The Kullback-Leibler (KL) divergence, D(f||g), between two distributions,</span>
0018 <span class="comment">% f(x) and g(x) is also known as the &quot;relative v_entropy&quot; of f relative to g.</span>
0019 <span class="comment">% It is defined as &lt;log(f(x)) - log(g(x))&gt; where &lt;y(x)&gt; denotes the</span>
0020 <span class="comment">% expectation with respect to f(x), i.e. &lt;y(x)&gt; = Integral(f(x) y(x) dx).</span>
0021 <span class="comment">% The KL divergence is always &gt;=0 and equals 0 if and only if f(x)=g(x)</span>
0022 <span class="comment">% almost everywhere. % It is not a distance because it is not symmetric</span>
0023 <span class="comment">% between f and g and also does not satisfy the triangle inequality.</span>
0024 <span class="comment">% It is normally appropriate for f(x) to be the &quot;true&quot; distribution and</span>
0025 <span class="comment">% g(x) to be an approximation to it. See [1].</span>
0026 <span class="comment">%</span>
0027 <span class="comment">% This routine calculates the &quot;variational approximation&quot; to the KL divergence</span>
0028 <span class="comment">% from [2] that is exact when f and g are single component gaussians and that is zero</span>
0029 <span class="comment">% if f=g. However, the approximation may be negative if f and g are different.</span>
0030 <span class="comment">%</span>
0031 <span class="comment">% Refs:</span>
0032 <span class="comment">% [1]    S. Kullback and R. Leibler. On information and sufficiency.</span>
0033 <span class="comment">%       Annals of Mathematical Statistics, 22 (1): 79�86, 1951. doi: 10.1214/aoms/1177729694.</span>
0034 <span class="comment">% [2]    J. R. Hershey and P. A. Olsen.</span>
0035 <span class="comment">%       Approximating the kullback leibler divergence between gaussian mixture models.</span>
0036 <span class="comment">%       In Proc. IEEE Intl Conf. Acoustics, Speech and Signal Processing, volume 4,</span>
0037 <span class="comment">%       pages IV�317�IV�320, Apr. 2007. doi: 10.1109/ICASSP.2007.366913.</span>
0038 
0039 <span class="comment">%       Copyright (C) Mike Brookes 2012</span>
0040 <span class="comment">%      Version: $Id: v_gaussmixk.m 10865 2018-09-21 17:22:45Z dmb $</span>
0041 <span class="comment">%</span>
0042 <span class="comment">%   VOICEBOX is a MATLAB toolbox for speech processing.</span>
0043 <span class="comment">%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html</span>
0044 <span class="comment">%</span>
0045 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0046 <span class="comment">%   This program is free software; you can redistribute it and/or modify</span>
0047 <span class="comment">%   it under the terms of the GNU General Public License as published by</span>
0048 <span class="comment">%   the Free Software Foundation; either version 2 of the License, or</span>
0049 <span class="comment">%   (at your option) any later version.</span>
0050 <span class="comment">%</span>
0051 <span class="comment">%   This program is distributed in the hope that it will be useful,</span>
0052 <span class="comment">%   but WITHOUT ANY WARRANTY; without even the implied warranty of</span>
0053 <span class="comment">%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the</span>
0054 <span class="comment">%   GNU General Public License for more details.</span>
0055 <span class="comment">%</span>
0056 <span class="comment">%   You can obtain a copy of the GNU General Public License from</span>
0057 <span class="comment">%   http://www.gnu.org/copyleft/gpl.html or by writing to</span>
0058 <span class="comment">%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.</span>
0059 <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span>
0060 [kf,p]=size(mf);
0061 <span class="keyword">if</span> isempty(wf)
0062     wf=repmat(1/kf,kf,1);
0063 <span class="keyword">end</span>
0064 <span class="keyword">if</span> isempty(vf)
0065     vf=ones(kf,p);
0066 <span class="keyword">end</span>
0067 fvf=ndims(vf)&gt;2 || size(vf,1)&gt;kf;       <span class="comment">% full covariance matrix vf is supplied</span>
0068 
0069 <span class="comment">% First calculate vf covariance matrix determinants and precision matrices</span>
0070 <span class="comment">% and then the individual KL divergences between the components of f</span>
0071 
0072 klff=zeros(kf,kf);      <span class="comment">% space for intra-a KL negative divergence</span>
0073 ixdf=1:kf+1:kf*kf;      <span class="comment">% indexes of diagonal elements in kf*kf matrix</span>
0074 ixdp=(1:p+1:p*p)';         <span class="comment">% indexes of diagonal elements in p*p matrix</span>
0075 wkf=ones(kf,1);
0076 <span class="keyword">if</span> fvf                  <span class="comment">% vf is a full matrix</span>
0077     dvf=zeros(kf,1);    <span class="comment">% space for log(det(vf))</span>
0078     <span class="keyword">for</span> i=1:kf
0079         dvf(i)=log(det(vf(:,:,i)));
0080     <span class="keyword">end</span>
0081     <span class="keyword">for</span> j=1:kf          <span class="comment">% calculate KL divergence between all mixtures and mixture j</span>
0082         pfj=inv(vf(:,:,j));
0083         mffj=mf-mf(j(wkf),:);
0084         pfjvf=reshape(pfj*reshape(vf,p,p*kf),p^2,kf);       <span class="comment">% pf(:,:,j)* all the vf matices</span>
0085         klff(:,j)=0.5*(dvf(j)-p-dvf+sum(pfjvf(ixdp,:),1)'+sum((mffj*pfj).*mffj,2));
0086     <span class="keyword">end</span>
0087 <span class="keyword">else</span>                <span class="comment">% vf is diagonal</span>
0088     dvf=log(prod(vf,2));
0089     pf=1./vf;
0090     mf2p=mf.^2*pf';
0091     mf2pd=mf2p(ixdf);       <span class="comment">% get diagonal elements</span>
0092     klff=0.5*(dvf(:,wkf)'-dvf(:,wkf)+vf*pf'-p+mf2p+mf2pd(wkf,:)-2*mf*(mf.*pf)');
0093 <span class="keyword">end</span>
0094 klff(ixdf)=0; <span class="comment">% force the diagonal elements to exact zero</span>
0095 <span class="keyword">if</span> nargin&lt;4
0096     d=0;
0097     klfg=klff;
0098 <span class="keyword">else</span>
0099     [kg,pg]=size(mg);
0100     <span class="keyword">if</span> pg~=p
0101         error(<span class="string">'GMMs must have the same data dimension (%d versus %d)'</span>,p,pg);
0102     <span class="keyword">end</span>
0103     <span class="keyword">if</span> nargin&lt;6 || isempty(wg)
0104         wg=repmat(1/kg,kg,1);
0105     <span class="keyword">end</span>
0106     <span class="keyword">if</span> nargin&lt;5 || isempty(vg)
0107         vg=ones(kg,p);
0108     <span class="keyword">end</span>
0109     fvb=ndims(vg)&gt;2 || size(vg,1)&gt;kg;       <span class="comment">% full covariance matrix vg is supplied</span>
0110 
0111     <span class="comment">% Calculate vg covariance matrix determinants and precision matrices</span>
0112     <span class="comment">% and then the individual inter-component KL divergences between components of f and g</span>
0113 
0114     klfg=zeros(kf,kg);      <span class="comment">% space for inter-a-b KL negative divergence</span>
0115     wkg=ones(kg,1);
0116     <span class="keyword">if</span> fvb                  <span class="comment">% vg is a full matrix</span>
0117         dvg=zeros(kg,1);    <span class="comment">% space for log(det(vg))</span>
0118         pg=zeros(p,p,kg);   <span class="comment">% space for inv(vg)</span>
0119         <span class="keyword">for</span> j=1:kg
0120             dvgj=log(det(vg(:,:,j)));
0121             dvg(j)=dvgj;
0122             pgj=inv(vg(:,:,j));
0123             pg(:,:,j)=pgj;
0124             mfgj=mf-mg(j(wkf),:);
0125             <span class="keyword">if</span> fvf              <span class="comment">% vf and vg are both full matrices</span>
0126                 pgjvf=reshape(pgj*reshape(vf,p,p*kf),p^2,kf); <span class="comment">% pg(:,:,j)* all the vf matices</span>
0127                 klfg(:,j)=0.5*(dvgj-p-dvf+sum(pgjvf(ixdp,:),1)'+sum((mfgj*pgj).*mfgj,2));
0128             <span class="keyword">else</span>                <span class="comment">% vf diagonal but vg is full</span>
0129                 klfg(:,j)=0.5*(dvgj-p-dvf+vf*pgj(ixdp)+sum((mfgj*pgj).*mfgj,2));
0130             <span class="keyword">end</span>
0131         <span class="keyword">end</span>
0132     <span class="keyword">else</span>                        <span class="comment">% vg is diagonal</span>
0133         dvg=log(prod(vg,2));    <span class="comment">% log(det(vg))</span>
0134         pg=1./vg;               <span class="comment">% precision matrix pg = inv(vg)</span>
0135         mg2p=sum(mg.^2.*pg,2)';
0136         <span class="keyword">if</span> fvf                  <span class="comment">% vf a full matrix, vg diagonal</span>
0137             vav=reshape(vf,p^2,kf);
0138             klfg=0.5*(dvg(:,wkf)'-dvf(:,wkg)+vav(ixdp,:)'*pg'-p+mf.^2*pg'+mg2p(wkf,:)-2*mf*(mg.*pg)');
0139         <span class="keyword">else</span>                    <span class="comment">% vf and vg both diagonal</span>
0140             klfg=0.5*(dvg(:,wkf)'-dvf(:,wkg)+vf*pg'-p+mf.^2*pg'+mg2p(wkf,:)-2*mf*(mg.*pg)');
0141         <span class="keyword">end</span>
0142     <span class="keyword">end</span>
0143     d=wf'*(<a href="v_logsum.html" class="code" title="function y=v_logsum(x,d,k)">v_logsum</a>(-klff,2,wf)-<a href="v_logsum.html" class="code" title="function y=v_logsum(x,d,k)">v_logsum</a>(-klfg,2,wg));
0144 <span class="keyword">end</span>
0145 
0146</pre></div>
<hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>